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ABSTRACT 

We investigate the damping of neutron star r-modes due to the presence of a vis- 

■ cous boundary (Ekman) layer at the interface between the crust and the core. Our 

study is motivated by the possibility that the gravitational-wave driven instability 
of the inertial r-modes may become active in rapidly spinning neutron stars, eg. in 

\ low-mass X-ray binaries, and the fact that a viscous Ekman layer at the core-crust 

■ interface provides an efficient damping mechanism for these oscillations. We review 

various approaches to the problem and carry out an analytic calculation of the effects 
due to the Ekman layer for a rigid crust. Our analytic estimates support previous 

j. | \ numerical results, and provide further insight into the intricacies of the problem. We 

f~>) ■ add to previous work by discussing the effect that compressibility and composition 

\Q \ stratification have on the boundary layer damping. We show that, while stratification 

■ is unimportant for the r-mode problem, composition suppresses the damping rate by 

about a factor of two (depending on the detailed equation of state). 

6 

£3 , 1 INTRODUCTION 

(Z3 ■ 

■ During the past few years we have witnessed a renewed interest in neutron star oscillations. One trigger for the recent 
I activity was the discovery that inertial modes in rotating fluid stars are generically unstable under the infl uence of gr avita- 

• ^ ■ tional radiation (a specific family of axial inertial modes, the so-called r-modes, being the most relevant) jAnderssonlll99 §1 ; 
r> \ iFriedman fc Morsinkll993l . Potential astrophysical implications of this instability have been extensively discussed in the lite r- 
ature, particularly in the context of gravitational wave observations (for a detailed review see I Andersson fc Kokkotasl l)200 lfl ). 
Anoth er motivation for work in this area is provid ed by the possible observations of crustal oscillations following magnetar 
flares dlsrael et al|l2005l : IWatts fc StrohmaveriEoOrl) . 

Soon after the discovery of the r-mode instability, it became clear that many additional pieces of physics need to be 
considered if one wants to draw "reliable" astrophysical conclusions. Since much of the required physics is poorly, or only 
partially, known this effort is ongoing. At the present time the most important damping mec hanisms that counteract the gr owth 
of an unstable mode are t hought to be i) viscous boundary layers a t the core-crust interface ijBildsten fc Ushomirskv l2000h and 
their magnetic analogue (Mendell 2001; Kinnev & Mendell 2003), ii) the bulk viscosity due to weak interactions in a hyperon 
core (iLindblom fc Owenl2002l : lNavvar fc Owe rl2005|), and iii) superfluid mutual friction due to s cattering of electrons off of the 
vortices in the neutron superfluid in the core J]VIendelllll99ll : lAndersson. Siderv fc Comerll2006l) . At the same time the mode 
ampl itude is expected to saturate due to nonlinear mode-mode coupling jArras^t all l2003l: iBriiik" Teukolskv fc Wassermanl 
2004). As the damping of unstable r-modes could be an issue of astrophysical relevance these mechanisms have attracted some 
attention. 

At the time of writing, the available results indicate that for mature neutron stars the key damping agent may be 
the friction associated with the crust interface. The crust forms when the star cools below 10 10 K or so, a few weeks to 
months after the star is born. Shortly after this the heavy constituents in the core (neutrons, protons, hyperons) are ex- 
pected to form superfluid /superconducting states. This leads to the suppression of, for e xample, the hyperon bulk viscosity 
ijflaensel. Levenfish fc Yakovlevll2002l : iReisene gger fc Bonacioj| 2f)03l : iNavvar fc Owen1l200fil) and the emergence of mutual fric- 
tion. The estimates for superfluid mut ual friction j M^n^j]^^^l|)indicate that, even though this mechan ism may suppress the 
insta bility in the f-modes completely l|Lindblom fc Mendell! 1 1993 ). it is not very efficient for r-modes ijLindblom fc Mendell! 



200c 

We therefore focus our attention on the role of the crust, which as a first approximation can be modelled as a rigid 
spherical "container" enclosing the neutron star fluid. Such configurations have been studied in several classic fluid dynamics 
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papers (see the standard textbook bv lGreensparJ Jl968l) for a detailed discussion of rotating fluids an d a wealth of references). 
In the context of r-mode damp ing the effect of a rigid crust was first d iscus sed some years ago by lfBildsten fc Us homirskvl 
i200d) . followed by the work of I Andersson et all feOQOh . iRieutorcj feOQlh and Lindblom. Owen fc Ushomirskvl feoool) . In this 
problem the core fluid oscillations are damped because of the formation of an Ekman layer at the base of the crust and 
the associated Ekman circulation. The result ing damp i ng is rapid enough to compete effectively with the gravitational- wave 
driving of the r-modes. Subsequent work bv iMendelll (|200lh indicates that when a magnetic field is present in the Ekman 
layer the damping timescale may become even shorter. Thus the dissipation due to the presence of the crust-core interface 
plays a key role for the r-mode instability. 

Our aim in this paper is to study the interaction between the oscillating fluid core and a rigid crust (using mainly 
anal ytical tools). A more r efined model, which considers an elastic crust, is discussed in an accompanying paper (see also 
iLevin fc Ushomirskvl J200lh ). Even though this problem has been discussed in some detail already, we have good motivation 
for revisiting it. It would be a mistake to think that the effect of the Ekman layer on an unstable r-mode is well understood. 
While we may know the answer for the simplest model problem of a viscous fluid in a (more or less) rigid container, a number 
of issues remain to be investigated for realistic neutron stars. The presence of multiple particle species and the a ssociated 
comp osition gradients may affect the Ekman circulation jAbnev fc Epstemlll996l) . the magnetic field boundary layer (IMendelll 
2001) depends strongly on the internal field structure (which is largely unknown) and whether the core forms a superconductor 
or not, once the core is superfluid one would need to account for multifluid dynamics as well as possible vortex pinning, and so 
on. To prepare the ground for an assault on the latter two problems, we want to obtain a better understanding of the nature 
of the viscous boundary layer and the methods used to study the induced damping of an oscillation mode. 

We will focus on the problem of a single r-mode in the core (a slight restriction, but the generalisation to other cases 
is not very difficult once the method is developed). In order for an analytic calculation to be tractable, a simple model for 
a neutron star is a prerequisite. Hence, we first consider a uniform density, slowly rotating, Newtonian star and reproduce 
familiar results for the Ekman problem. In doing this, we develop a new scheme which should prove useful in more complicated 
problem settings. Our study also serves as an introduction to boundary layer techniques for the reader who is not already 
familiar with the tricks of the trade. We then compare and contrast the different methods that have been used to investigate 
the Ekman layer problem. This discussion serves to clarify the current understanding of the in duced r-mode damping r ate. 
Finally, we extend our analysis to compressible fluids with internal stratification. The results of lAbnev fc Epstein! Jl99fill for 
the spin-up problem suggest that both these features will have a strong effect on the Ekman layer damping. However, we prove 
this expectation to be wrong. For a toroidal mode, like the r-mode, we find that stratification is irrelevant. T he compressibility, 
on the other hand, does play a role but we demonstrate that the effect is much weaker than indicated bv lAbnev fc Epstein! 
( 1996() . Rather than leading to an exponential suppression, the compressibility is estimated to reduce the r-mode damping by 
about a factor of two. 



2 FORMULATING THE EKMAN PROBLEM IN A SPHERE 

As a first application of our boundary layer scheme we consider the E kman layer prob lem in a spherical container. This 
is an instructive exercise as the solution to this problem is well known jGreenspanlll968f) . The so-called "spin- up problem" 
forms a part of classic fluid dynamics and was s tudied in a numb e r of papers in the 1960s and 70s. We found the work by 
iGreenspan fc Howard! <|l96,^ : loreenspanl (|l964h : Iciark et all l|l97lh : lAcheson fc Hvdel <|l97:t : iBenton fc Clark! l)l974 partic- 
ularly useful. The Ekman problem is also relevant for a discussion of unstable r-modes. In fact, it is reasonable to expect that 
the associated r-mode damping timescale is of the same order of magnitude as the spin-up timescale in a rigid sphere since 
an inertial mode has frequency proportional to the rotation rate of the star, a tx Q. 

In building a simple model of a neutron star with a rigid crust, we assume a core consisting of an incompressible, uniform 
density, viscous fluid surrounded by a rigid boundary. Furthermore, we work in a slow-rotation framework (which greatly 
simplifies the mathematical analysis as the star is retains its spherical shape at O(0,)). Finally, we make the so-called Cowling 
approximation, i.e. neglect any perturbation in the gravitational potential. These simplifications are obviously not necessary 
in a numerical analysis of the problem, but (in our view) the benefits of having an analytic solution far outweighs the value of 
a numerical solution to the complete problem. One can argue in favour of some of our assumptions, since most astrophysical 
neutron stars are slowly rotating (in the sense that their spin is far below the break-up rate) and the inertial modes which 
we are interested in are such that the density (and hence gravitational potential) perturbations are higher order corrections. 

Working in the rotating frame, our dynamical equations are the (linearised) Euler and continuity equations, 

d t v + 2fi x v = -Vh + iA7 2 v (1) 
V • v = (2) 



The enthalpy h, which is defined as 
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Vh = -Vp (3) 
P 

has been introduced for later convenience. In addition to the equations of motion, we have two boundary conditions for the 
core flow: (i) regularity at r — and (ii) a "no-slip" condition v(_R c ) = at the crust-core interface r = R c . The role of 
viscosity is to ensure that the second condition is satisfied by modifying the inviscid solution in a thin boundary layer near 
the crust. This is known as the Ekman layer, and it has characteristic width Se ~ \J~ujVL ~ 1 cm for typical neutron star 
parameters. Our aim in the following is to find a consistent solution to the problem when the solution in the bulk of the 
interior (the inviscid flow) corresponds to a single I — m r-mode. 

Sc hematically we assume that the solution takes the following form (a standard uniform expansion in boundary layer 
theory lender OrzaJlimltt: 

V = V +V +l^ 1/2 [v 1 +V 1 ] +0(v) (4) 
and similarly for the perturbed enthalpy Sh and the mode frequency a (we assume that all perturbations behave as exp(icrt)) 
Sh = 5h° + 5h° + v 1/2 [Sh 1 + 5h 1 ] + 0{v) (5) 

a = (T + v 1/2 is + Oiy) (6) 

The complex-valued frequency correction s is of particular importance, since its real part corresponds to the Ekman layer 
dissipation rate. The idea behind the expansion (0J is that the inviscid solution v° is modified by v° in the boundary layer in 
such a way that the no-slip condition is satisfied. This then leads to a change in the interior motion, represented by v 1 , which 
requires another boundary layer correction v 1 and so on. In principle, this leads to an infinite hierarchy of coupled problems. 
In this picture, the various boundary /crust-induced corrections to the basic inviscid flow fall into two categories, namely, 
boundary-layer corrections (always denoted by tildes in this paper) which are characterised by rapid (radial) variation [here 
d r ~ 0(^ _1//2 )] and secondary "background" corrections which vary smoothly throughout the star. Based on this different 
behaviour, when Q is inserted in Jy and each of these equations splits into two components; one containing only 
boundary layer quantities and another comprising the background quantities. These two families of equations are not entirely 
independent, since the boundary conditions communicate information between boundary layer and background quantities. 
Assuming that the solution to the inviscid problem takes the form of an axial r-mode we get 

d t v° -2ncos8v° = --dgSh (7) 
r 

d t vi + 2fJcos0vg = 3_A<5h° (8) 

20sin(9u° = d r 5h° (9) 
We also know that the I = m r-mode solution corresponds to a frequency 

and can be expressed in terms of a stream function U° as 

v o = ^-id v U° + ^idoU (11) 

r sin ff r 

where the phase is chose for later convenience, and 

( r \ m+l 
RJ (12) 
(normalised to unity at the base of the crust). 

The corresponding boundary layer correction is a solution to 

d t vg -2ficos0«° = udlvg (13) 

dtthp + 2Q cos 8vg — ud 2 v^ (14) 

2Q.sin6v° p = drSh 1 (15) 

d r Sh° = (16) 

We also need to satisfy the continuity equation which leads to uj? = 0, and 

v l/2 rd r v l r + -±-de(smev° 9 ) + -^d v vl = (17) 
smt/ suit/ 

Note that we have only kept the highest radial derivative of the various variables in the boundary layer, and that we have also 
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neglected angular derivatives compared to radial ones. Finally, the solution to the v° problem must be such that the no-slip 
boundary condition is satisfied, i.e. we must have 

v°+v°=0, at r = R c . (18) 

At the next level, we determine the change in the interior motion induced by the presence of the Ekman layer. This means 
that we must solve, 

dtvl — 20 cos Qv\, — — -dgSh 1 + svg (19) 
r 

dtvl + 20 cos Ov] + 20, sin Bv\ = ^d^Sh 1 + svl (20) 

d t vl — 20 sin 6v]p — —drSh 1 (21) 
together with the continuity equation 

V ■ v 1 = (22) 
and a boundary condition 

Vr + Vr=0 (23) 

A solution to this problem should provide the frequency correction s, which encodes the rate at which the presence of the 
Ekman layer alters the core flow. That is, the induced damping of the r-mode which we are interested in. 



3 SOLUTION IN TERMS OF A SPHERICAL HARMONICS EXPANSION 

So fa r, our calculation follows closely the standard treatment of the spin-up problem jGreenspan fc HowardllT963t iGreenspanI 
1968). At this point we deviate, by expanding the velocity components in spherical harmonics. We have several reasons for 
doing so. Firstly, the standard approach (which will be reviewed in Section|lJ leads to a boundary layer flow that is explicitly 
singular in particular angular directions. Although one can demonstrate that this is not a great problem, and that the predicted 
Ekman layer dissipation can be trusted, it is not an attractive feature of the solution. We want to use an approach that avoids 
such singular behaviour. Secondly, the standard approach is somewhat limited in that it may not easily allow us to consider 
a differentially rotating background flow (or an interior magnetic field with complex multipolar structure). Our approach will 
allow such generalisations, should they be required. 

3.1 The boundary layer corrections 

We first concentrate on the correction to the fluid velocity in the Ekman layer. We need to solve 11311 - 114(1 and then determine 
v r from 1171 . If we are mainly interested in determining the damping timescale due to the presence of the Ekman layer, the 
induced radial flow will play the main role. 
Assuming that v° can be written 

v = — e r + - dgV t—t ) eg + - + idgll e v (24) 



r r V smf / r V sm 

i.e. as a general decomposition into polar (W, V) and axial (U) perturbations, and that our solution depends on ip as e tmip , 
we get from H'M) 

m cos 6V + sin 8 cos 6dgU = -j^o + ivti?)(sw.0daV + mil) (25) 
Similarly, equation (11411 leads to 

sin6»cos6>d 9 V' + mcos6>[/ = — (<r + ivdr)(mV + sm6d e U) (26) 

We can combine these two equations to get 
- — (ao +iud*)VeU = mil - cos8X7 2 g V + sinOdgV (27) 
r^(<ro + wdr)VeV = -mV + cos 9X7 2 e U -sin 6dgU (28) 
where 
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is the angular part of the Laplacian on the unit sphere. 

If we now expand the various functions in spherical harmonics, i.e. U = ^ ( UiY™ etc, and then use their orthogonality 
together with the standard recurrence relations, 

coBWr = Qi+rF^i + QiYi-x (30) 
smOdeYT = IQ1+1YR1 - (I + l)QiY£i (31) 
where 



(I + m)(l — m) 

Ql = V(2i-1)(2/ + 1) (32) 
we arrive at the two coupled equations 
n(n + 1) 



2<> 



(<to + iud r )U n — mU„ — (n — l)(n + l)Q n V n -i — n(n + 2)Q n +iV n +i = (33) 



+ 20 + (g ° + iua %)V n +i - mV„ + i - n(n + 2)Q n+1 U n - (n + l)(n + 3)Q w+2 £/n+2 = (34) 

Since we know that we must match our solution to the inviscid r-mode solution at r — R c , we now assume that all 
toroidal components apart from U m vanish. Strictly speaking, this assumption is inconsistent. This can be easily seen by 
setting n — m^ra + 2in the first of equations above. In principle, o ne would expect an infinite sum of contributions. This 
problem is similar to that for an inertial mode in a compressible star l|Lockitch fc Friedmanl lT999 ) . From a practical point of 
view, however, one must truncate the sum at some level in order to close the system of equations. Just like in the inertial-mode 
problem, one would expect the solution to converge as further multipole contributions are accounted for. This then provides 
a simple way to estimate the error induced by the truncation. Adopting this strategy, we will later compare our results to 
ones obtained by including more axial terms (U m +2, Um+i and so on) in the boundary layer flow. 

After setting n = m we get (since Q m = 0) 

m{m + 1) 

(m + l)(m + 2) 
20, 

Using the standard ansatz [C/ m , Vm+l] oc e Ar for a differential equation with constant coefficients, we can show that a 
solution will exist if 

2,4 4ifi ^ 2 4m(m + 2)Q 2 Q 2 m+1 

VA H -, ' TTo = (37) 



■((Jo + iud r )Um — mUm — m(m + 2)Q m+ iV m +i = (35) 
(cr + iudr)V m +i - mVm+i - m(m + 2)Q m+1 U m = (36) 



(m+l)(m + 2) (m + 1) 

i.e. for the two roots 
>2 2iQ 



(m + l)(m + 2) 
where we have defined 



(l±a) (38) 



a = y/m(m + 2)3Q* m+1 + l (39) 
The solutions correspond to eigenvectors 

V m+ i = - , m 2i*n Um (40) 

(m + 2)^Q m+ i 

In order to satisfy the prescribed boundary conditions at the crust-core interface, we must have a vanishing velocity at 
r — R c . Writing the solution as 

U m = Ae Al(r - flc) + Be X ^ r - Ra) (41) 
where 



Ai, 3 = (l±t)^ p (42) 
" v(m + l)(m + 2) 

we must have 

Um (Rc) = U° m (R c ) — ► A + B = 1 (43) 
Furthermore, we also require that 

V m+1 (R C ) = , i 2n [(1 + a) A + (1 - a)B] = (44) 
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which leads to 

A = -— and B = — (45) 

Finally, we can use this solution together with (11711 to determine v T . Assuming that v 1 also takes the form (12411 . we readily 
find that 

v 1/2 d r {rW 1 t+1 ) = (m+ l)(m + 2)VWi (46) 
which can be integrated to give Hence, we have obtained all leading-order boundary layer variables. 

3.2 The induced interior flow 

Having determined the appropriate solution in the Ekman layer, we now want to calculate the induced flow in the interior 
and deduce the damping rate of the r-mode. By combining (11911 and 121 It . and again using a solution of form 1241 [albeit now 
without the tildes] for v 1 , we can show that we must have 

( ao V^t/ + 2mQU) - 2ttcos0S7gV + 2ftsmed e V - 2tt{sm0d e W + 2cos6>W] = -is\7 2 g U° (47) 

Expanding once again in spherical harmonics and using the fact that the right-hand side vanishes unless I = m, we immediately 
arrive at 

2(m + 2)nQ m+1 V m+ i + 2ftQ m+1 W m+1 = is(m + 1)17° (48) 
and a solution 

ir s(m + l) o 1 „, 

Vm+1 = OCX i J U ™ -r^WWl 49 

2Q(m + 2)Q m +i m + 2 

Given this result, we can use the continuity equation to show that 

d r (rW m+1 ) = (m + l)(m + 2)V m +i = ^±12! ^ - (m+ l)HWi (50) 

(Note that the integration constant must vanish since we would otherwise have a solution that diverges at the origin.) 
The final step in our analysis corresponds to combining 1461 and 1501 in such a way that 

W m+1 (R C ) + W m+1 (R C ) =0 (51) 

By integrating the two solutions we find that 



1/2 r> Txr m + 1 

/ ' -RcWm + l = 



1 1 

Ai A2 



(52) 



m + 2 2Q m+ i a 
and 

W m+1 (Rc)= im + 1) ^ SQm+1 (53) 
After combining these and reshuffling, we have 



n (2m + 3) i(a 2 - 1) 



1 1 

Ai A2 



(54) 



Im s ; 



v 1 / 2 R c (m + 2)(m + l) a 
From this last equation we can extract the damping rate, 

Re s = ( " ) 2m + 3 ^TT+^T /T7 (55) 

yR 2 J v / 4(m + l)(m + 2) a 

We also find that the frequency correction due to the presence of the Ekman layer is 

(") 1/2 , 2m + 3 ^-^V^-l (56) 

\RcJ > /4(m + 2)(m + l) Q 

These expressions constitute the main new result of the first part of this paper. They have been derived by combining 
boundary-layer theory techniques with a decomposition in spheri cal harmonics. Pu tting nu mbers in, we find for m = 1 
that s w 2.92 — 0.602i. This should be compared to the results of lOreenspanl l|l968h (see also iRieutordl ^OOlh ). who finds 
s ss 2.62 — 0.259i. The quantity which we are mainly interested in, the real part of s which gives the damping rate of the core 
flow, differs from Greenspan's result by about 11%. This is not too bad given that we truncated the boundary layer solution 
at the level of a single multipole contribution. If we include further terms, we need to solve a system of equations larger than 
the "2 x 2" problem given by 1331 - 1341 . Including one more axial term (Um+2) in the boundary layer we end up with a "3 x 3" 
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truncation 


This work (Section 151 


Greenspan ("Section |41 


m = 2 2x2 


4.097 - 0.473i 


3.518 - 0.176i 


3x3 


3.451 + 0.376i 




4x4 


3.714 - 0.292i 




m = 3 2x2 
3x3 
4x4 


5.067 - 0.391i 
3.922 + 0.592i 
4.405 + 0.209i 


4.288 - 0.136i 



Table 1. Ekman dissipation rate s computed in the present work, compared with the standard approach in the literature iGreenspanl 
1968). 



problem and so on. As we include more multipoles in the description our damping rate approaches Greenspan's result, cf. the 
results listed in Tabled We see that if we consider the "4 x 4" problem we get results that agree with the standard ones, see 
the following section, to within 5% or so. This is certainly much more accurate than our knowledge of the various neutron star 
parameters. It should be noted that, by solving the problem in a "non-standard" way we have learned an interesting lesson: 
In order to represent the boundary layer flow accurately, we need to include a number of multipoles even though we have a 
very simple core flow. 



4 THE STANDARD SOLUTION TO THE EKMAN PROBLEM 

For completeness, and comparison with our calcul ation, we will now discuss the "traditional" solution of the Ekman prob lem. 
This solution is described in detail in the paper bv lOreensoan fc Howard! l)l963h as well as the book bv lOreensrjanl l)l96sF) and 
is, essentially, a boundary layer theory calculation without expanding in spherical harmonics. The formulation of the problem 
is identical to that described above. 

By combining I13H and 11411 we can easily show that we must have 

\y 2 dr — 2ioovdl — <7q + 4fi 2 cos 2 #] vg = (57) 

Despite really being a partial differential equation, this is a now viewed as a "constant coefficient" ordinary differential equation 
for vg as a function of r. The ansatz vg = e A '' R °' then leads to the indicial equation 

u 2 \ A - 2ia Q v\ 2 + 4fi 2 cos 2 9 - a 2 = (58) 
which has solutions 

A = ±^ 1/2 [i(cr ±2ficos6i)] 1/2 (59) 

In order for the boundary layer corrections to the velocity field to die away from the surface we must have Re A > 0. One can 
show that the two solutions that satisfy this criterion correspond to 

1/2 " 



Al 2 — 



m+1 



± cos 6 



! l/(m + l)±co8g . 
+ l/(m + l)±cos6»| I 



(60) 



where we have i) used the I = m r-mode frequency, and ii) introduced A = •jQ./vX. 
Given a solution for vg we find from 1 1 H that 

± iv v = vg (61) 

where the upper (lower) sign corresponds to Ai (A2). The leading order correction to the velocity field due to the presence of 
the Ekman layer can now be written 

vg = Ae Xl{r ' R ' ] + Be X2{r - R ' ] (62) 
v v = -iAe Xl(r - Ra) +iBe X2{r - Rc) (63) 

The coefficients A, B here should not be confused with those appearing in Sectional (although in both cases they have similar 
roles). This solution needs to be matched to the inviscid solution in such a way that the no-slip condition at r — R c is satisfied. 
Matching to the I — m inviscid r-mode leads to 

A = 2^ke Y - ~ 2k deY - (64) 
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This completes the solution in the boundary layer. We want to infer the induced radial flow in order to estimate the Ekman 
layer damping of the core r-mode oscillation. The required relation is best expressed in terms of a "stretched" boundary layer 
variable 

C = yf(iic-r) (66) 

Then the 0(1) continuity equation leads to 
-i 1 



sinfl sinS 



or 



V^Rl 



-^—rdg(sm9v e ) + -^—p:d v v v 
smO smO 

,1/2 



ld( 



This is to be matched to the 0(y ' ) core solution in such a way that 
lim v r = lim v r 

r — >R C Q — >oo 

Working things out for our boundary layer solution we find that 



lim vl — , , 



1 



im(rn + 1) 
2R C 



1 

17 



i 

T 2 



Y" 



A 



A? 



B 



d e \ 2 

A 2 



Using the induced core solution from before, i.e. 1)53^ . we have 



(m+ l) 2 s 



-Y, 



m+l 



r ™ c Vr ~ 2{2m + 3)R c QQ m+1 
and we can determine the Ekman layer dissipation rate for the r-mode from 



(m + l) 2 s 



-Y 



2(2m + 3)Q m+ r m+1 V /»'r I 2 
The final step is to use orthogonality of the spherical harmonics, i.e. evaluate s from 
2(2m + 3)Q„ 




Q, J im(m + 1) 



i 

17 



i 



Y" 



AR, 



BR C 



(m + l) 2 



where 



3 = 2tt 



sin 6Y2 



im(m + 1) 



1 

17 



i 

17 



± IT, 



AR, 



A? 



BR C 



d g \2 

A 2 



dO 



(67) 



(68) 



(69) 



(70) 



(71) 



(72) 



(73) 



(74) 



The required integral can be much simplified if we use the symmetry of the integrand, see the discussion in the Appendix or 
alternatively iLiao. Zhang fc Earnshawl teOQll) . Thus we get 



J = 2-Ki 



sin 6 



ym yr, 
! 1 m+l 1 rr, 



-cMYT+iKr) + (dgYZ +1 )(deYZ) } d9 



(75) 



The obtained values for the damping rate s are given in Tabled 

The solution to the Ekman problem described in this Section is straightforward but somewhat limited, in the sense that 
it relies on the simplicity of equation 1571 . In a more general context, this approach may no longer be possible. Moreover, as 
we already pointed out, this solution exhibits a singular behaviour for certain angles (when A = 0). From an intuitive point 
of view one would clearly not expect a physical quantity, in this case the velocity field of the fluid, to exhibit singularities in 
certain angular directions. Hence, it is useful to discuss these singularities in more detail. To unveil the origin of the singular 
solution we need to return to Eq. Q. Key to the boundary layer approach was the assumption that the viscous term on the 
right-hand side is relevant only in a narrow region in the vicinity of (in our case) the neutron star crust. The idea is that 
the inviscid solution v° is modified by a correction v° in the boundary layer in such a way that a no-slip condition can be 
imposed. The boundary-layer corrections are characterised by rapid (radial) variation [d r ~ 0(i/ -1 ' 2 )]. Thus, neglecting the 
angular derivatives in the Laplacian in the right-hand side of the Navier-Stokes equation we can derive Eqs 11311 - 114(1 . If we 
then conveniently "forget" that we are actually solving partial differential equations and treat 6 as a "constant" then we 
arrive at 157^ . At the end of the day, the solutions we obtain are manifestly singular for particular values of 6. The problem 
stems from the assumption that the radial derivatives dominate the angular derivatives in the boundary layer. Without this 
assumption we would have to solve the original partial differential equation, the standard approach to which is (of course) 
separation of variables via an expansion in spherical harmonics (as in our method in Section However, given the singular 
solution derived above it is very easy to show that 
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Figure 1. We compare the solution of the 2x2 problem (dashed) to the standard "exact" solution (solid line). The data corresponds 
to parameters m = 2, E = v/2QR 2 = 10 -6 , ip = ir/3 and a typical angle 9 = 7r/4. The two solutions are clearly in good agreement. 
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Figure 2. We compare the solution of the 2x2 problem (dashed) to the standard "exact" solution (solid line). The data corresponds 
to parameters m = 2, E = vj1Q,B? = 10 — 6 , ip = n/3 and a distance away from the boundary given by x = r/R = 1 — 0.8V1?. The 
particular values of 9 for which the standard solution is singular are apparent in the figure. However, away from these angles we find the 
two solutions to agree well. 



dv e dv e , 

« -w (76) 

near the singular angular directions. Hence, this solution violates the assumptions made in the derivation which renders the 
calculation inconsistent. This argument shows that this solution is u nph ysical in the vicinity of these angular directions. 
Similar conclusions have been drawn bv lRieutord fc Valdettarcl Jl99?l) and Irlollerbach fc KerswelH Jl995l) . They consider the 
complete problem numerically and show that shear layers form in these specific angular directions. They also prove that the 
presence of these shear layers have a small effect on the damping rate. This is the key reason why both our method from 
Section El and the "exact" singular solution derived above, where the detailed shear layers are neglected, lead to consistent 
and reliable estimates for the damping due to the Ekman layer. A comparison of the two velocity fields, cf. Figure Q and 
[5] shows that the agreement is quite good away from the singular directions. Hence, it is not very surprising that the two 
calculations lead to similar estimates for the r-mode damping. It is worth noting that iHollerbach fc Kerswelll lll995t) argues 
that the relatively slow convergence of the correction to the oscillation frequency, cf. Tabled is due to the presence of the 
shear layers. 



5 ANOTHER SOLUTION: SURFACE INTEGRAL APPROACH 

There is yet another alternative for solvin gj the Ekman pro blem, based on direct integration of the Euler equat ion. This 
appro ach i s discussed in Greenspan's book <IGreensDan1ll968i) . With slight modifications, this method was used bv iRieutordl 
feOQll) and lLiao. Zhang fc Earnshawl J200lj) in their studies of inertial modes. 
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Returning to equation I12H , which determines the induced core flow, and reinstating the pressure and density instead of 
the enthalpy, we have 

ictov 1 + 20 x v 1 + -X/Sp 1 = sv° (77) 
P 

Multiplying by the complex conjugate v° (the conjugate is denoted by a bar), and using the (conjugate of the) inviscid 
equation 

ia v° + 20 x v° = -V8p° (78) 
P 

as well as the continuity equation V • v° = 0, we obtain, 

V ■ [SpV + 5pV] = ps|v°| 2 (79) 
Integration of this relation over the volume (V) of the star and use of the divergence theorem gives, 

/ r • [(JpV + [Sp^v 1 ] dS = s I p\v°\ 2 dV (80) 
J s Jv 

where S represents the surface at the core-crust interface. We know that the inviscid r-mode problem is such that 

f-v° = (81) 



which means that we arrive at the final expression 

(82) 



J>|v°|W 

Now we can combine 1821 with the matching condition for the induced flow 1691 1. This provides a connection to the 
calculation in the previous section. We can easily determine 8p° for the inviscid flow. From the r-mode solution to the velocity 
field we find that 

n 2m t r \ " n + 1 

Sp° = -i — pO — Q m+1 Y m+1 (83) 

m + 1 \RcJ 

To determine the damping rate we also need the volume integral (essentially the mode energy) 
Putting all these ingredients together we find that (using definitions from the previous section) 



2( 2m + 3)Q m+ i []!_ 
(m + 1) 2 V Ri 

which is in agreement with the result, Eq. 1731 . of the previous section. 

This final approach to the problem is obviously very elegant, and much less involved that the previous two calculations 
we have presented. Unfortunately, the surf ace integral approach r elies upon the assumption of incompressible flow. This 
assumption will not hold for stratified stars jAbnev fc Epstem]ll99cj) . see Section |7| or indeed stars with an internal magnetic 
field. Hence, we do not expect this method to be of much use for more complex (realistic) problem settings. 



6 COMPARING DIFFERENT RESULTS 

Before we move on to consider the more complex problem of a compressible star, it is meaningful to compare and contrast the 
different results for the Ekman layer damping of r-modes in the current literature. Let us do thi s in chronological order. That th e 
viscous boundary layer would provide an efficient damping mechanism was first suggested bv lBildsten fc Ushomirskvl J2000l) 
They estimated the associated damping rate by using the familiar plane-parallel oscillating plate result from lLlmdlaTfcljrfsnha 
(1987]). Their result is not too far from the more detailed results derived in this paper, even though it does not account for 
the geometrical factors due to the sphericity of the star and the nature of the r-mode fluid flow. This is not very surprising 
since the geometric factors are of order unity as long as the thickness of the boundary layer is considerably smaller than the 
size of the star. At a local level, the bound ary layer flow is quite similar to that in the plane-parallel case. 

In his comment on these first estimates. lRieutordl l)200ll) extracted l|73|l from Green span's work. As we have seen, this result 
has the relevant geometric factors etcetera, but it is based on a uniform density model. iLindblom. Owen fc Ushomirskvl l|2000fl 
extended the standard method discussed in Section [1] to allow for non-uniform density models, and derived the boundary 
layer-modified r-mode flow. Then, using the viscous shear tensor they arrive at the final result for the r-mode damping 
timescale, 
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t _ 2E _ 2 m+3/2 (m + l)! /2Q,Rip c f R ° p ( r \ 2m + 2 dr 



E 2mfi(2m+ l)!!X m \j t] c J p 



/ r \ -zm+z J r 



where rj c and p c are the viscosity coefficient and the density at the core-crust interface, respectively, and T m is an angular 
integral analogous to our J . The integral on the right-hand side is essentially the kinetic energy of the r-mode. In order to 
compare this expression to th e result quoted b y Rieutord, we should take p — p c = constant and work out the radial integral. 
This readily yields Eq. (11) in iRieutordl (|200lt ). Hence, the two results are consistent, as they should be. To compare with our 
result, Eq. (I73H . we need to use the relation between the two angular integrals. One can show that 

Raj= m( y a 3) " Q m+ A, (sr) 

Given this relation our result for the damping timescale 



(m+1)! pj£_}_ 

E 2(2m + 3)Q m+ i V VcV Re J { ' 

is also identical to Rieutord's result. From this comparison it is clear that any difference between the various detailed studies 
of the Ekman layer problem in the r-mode instability literature is due to the chosen stellar model. 

To quantify the dependence of the damping rate on (say) the mass distribution in the star is obviously an important task. 
ILindblom. Owen fc Ushomirskvl J2OOO) did this by performing calculations for a number of realistic supranuclear equations 
of state. Here we will choose a more pragmatic approach, which we think provides some useful insights. Let us first return to 
1861 1. but instead of taking p — p c we will, somewhat inconsistently, allow the two (constant) densities to be different. Then 
we see that, in the case of the m = 2 mode, 



/= r^m 6 £=^3.4ixio- 

J pc \RcJ Rc 7pc 



(89) 



R?p c 

Now compare this result to a star described by an n = 1 polytrope for which 

If we t ake p c = 1 .5 x 10 14 g/cm 3 and R c /R ~ 0.85, which is a fairly typical value according to Table 1 in lLindblom. Owen fc Ushomirskvl 
(2000), we find that 



/« 3.49 x KT 2 --— (91) 
R 6 p c 

From thi s comparison we see that one can quite easily bridge the difference between the constant density calculation and the 
results of lLindb lom 1 Owen fc Ushomirskyj <|2000 f). A s a further exercise, one can work out the integral by combining the data 
in Table iof lLhKibl'o'rn^o'wCT^ for M, R and R c with the n — 1 polytropic density profile. This leads to 

underestimates of the damping timescale results obtained by Lindblom et al for various equations of state by 10-30%. This 
is due to the fact that realistic equations of state are slightly softer than the n = 1 polytrope. Nevertheless, it shows that 
the polytropic model captures the main features of the problem. One should keep in mind that the use of realistic equations 
of state for Newtonian neutron star models is dubious. For a given central density (say) the mass and radius often differ 
considerably from the corresponding results in general relativity. This means that it may not be entirely consistent to use 
realistic equations of state in the present discussion. On the other hand, it is important to get an understanding of how the 
results change with the stellar parameters, and over what range the damping timescale for an unstable r-mode may vary. We 
think that the above discussion led to some useful comparisons that provide useful steps towards this understanding. Possibly 
the most important conclusion is that all current calculations are consistent. 



7 INCLUDING STRATIFICATION AND COMPRESSIBILITY 

So far, we have discussed the Ekman layer problem for incompressible fluids. Although this is not an unreasonable first 
approximation for neutron stars, it is well documented that both compress ibility and stratification d ue to either temperature 
or composition gradients may be important. Considering also the results of lAbnev fc Epstein! l|l99fl) for the spin-up problem, 
which suggest that the Ekman layer timescale is significantly altered by both stratification and compressibility, it is clear that 
we need to relax our assumptions. 
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7.1 Formulating the compressible problem 

We want to avoid the assumption of uniform density and include stratification and compressibility (while maintaining the 
calculation at O(0.) which means that the star is spherical). In a perturbed compressible, stratified star the equation of state 
is usually taken to have the form 

Ap = ^-Ap (92) 
P 

where A represents a Lagrangian variation and T is the adiabatic index. Written in terms of Eulerian perturbations (S) this 
means that, for our inviscid leading order flow we have 

!£ = g.-A.e (93) 
p Fp 

where the Schwarzschild discriminant is 

A s = P --f (94) 
P Fp 

(the pri mes represent radial derivatives). The Schwarzschild discriminant encodes information regarding composition gradients 
etcetera iReisenegger fc Goldreichlll992t) 

The question now is, how does this more complex equation of state affect the boundary layer calculation? At first sight 
one might think that the answer will be quite messy. Fortunately, this turns out not to be the case. Let us first consider the 
inviscid flow discussed in Section [5] Then we see that the introduction of the enthalpy was, indeed, a shrewd move. With the 
equations written in this way we do not need to worry about the relation between the perturbed density and pressure unless 
these variables are explicitly required. In our case, they are not. Thus we swiftly prove that the inviscid r-mode flow remains 
unchanged for a compressible model. 

The leading order boundary layer flow also requires little thought. As before we find v° = 5fP — 0. Given this, the 
remaining equations at this level are identical to those from the uniform density calculation. Hence Vg and 0° are obtained 
as in Section Finally, from the equation of state we find, 

£-f£-o <«> 

P Fp 

The complications associated with stratification and compressibility enter the problem at the level of the induced core 
flow. Expressed in terms of the enthalpy the Euler equations 119H - 12HI remain unaltered, but the continuity equation is now 
given by 

■ ^ P9 i , V7 i , is ^ s n ma\ 
too — — v r + V ■ v H v r = (96) 

1 P 1 P (TO 

where we have prefered to use the perturbed pressure in order to emphasize the presence of the stratification. Inspection of 
the Euler equations, together with the fact that s ~ 0(f2 1//2 ), implies that Vg „ ~ 0(f2 -1//2 ) and Sp 1 ~ O{0}^ 2 ). Hence, in the 
continuity equation the second and third term are 0(f2 _1,/2 ) while the first term is of order C(f2 3 ^ 2 ). At the same time, the 
fact that we are considering an r-mode flow means that the last term is also 0(f2 3//2 ). As we are working in the slow-rotation 
limit we keep only the former terms, and we get 

V ■ v 1 - Pi-vl = (97) 
Fp 

where g — p' / p. It is worth remarking that lAbnev fc Epstein! (^.996) make similar assumptions in their analysis of the spin-up 
problem. At this point we can infer that stratification has no effect (at leading order) on the Ekman damping rate of a 
r-mode. This is in contrast to the situation for a general inertial mode (and the spin-up problem) for which the radial velocity 
component is of the same order as the angular ones. Then the last term in I96H must be retained and stratification could be 
important. For our problem this turns out not to be the case, which may be a bit of a surprise. 

Our next task is to calculate the induced flow v 1 . Asssuming the usual decomposition of the velocity field we retain 1491 . 

i.e. 

t/1 _ is{m+l) o 1 w i rqs x 
At the same time the continuity equation gives, 

rdrWUi + [l - Wm+l = (m + l)(m + 2)V^ +1 (99) 
where we have defined the dimensionless "compressibility" , 
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Combining these equations we arrive at 



rdrWL. 



+ 



m + 2 —K 

R . 



wL 



is(m + if 
2QQ m +i 



'-- n 



(101) 



Once we solve this equation, we can replace 1531 in the analysis of the incompressible case and infer the damping rate of the 
r-mode. The problem is that, whereas the k = case is easily dealt with, the radial variation of k in a realistic model will 
typically require a numerical solution. 



7.2 Case study: The n = 1 polytrope 

To get an idea of how important the compressibility corrections may be we will consider the simple m odel problem of an 
n = 1 polytrope. Before doing this, let us make a few comments on the work of lAbnev fc Epstein! il996i) . In their analysis of 
the spin- up problem they arrive at an equation corresponding to IIUIH . They then make the assumption that k can be taken 
as constant. Since k ^> 1 in the region of the boundary layer the calculation simplifies considerably, and one easily shows that 
the compressibility leads to an exponential suppression of the Ekman layer damping. Working out the relevant algebra, we 
find a revised Ekman damping rate s given by 

(p \ 2m+3 
^) e-"-'* (102) 



s (2m + 3)! 

For canonical values m = 2, R c = 0.9-R and k — 10 (see below) we would have 

-wO.l (103) 

s 

T his suggests that the Ekm an layer damping is strongly suppressed in a compressible star. In fact, if we take k « 100 as 
in lAbnev fc Epstein! Jl996ll we would have a truly astonishing result where the boundary layer dissipation was completely 
irrelevant. The analysis is, however, wrong. 

To see what the problem is, we can work out n for the polytrope 19UII . This leads to 

k = — (1 — 7r cot -nx) w — — — (104) 
x 1 — x 

where x = r/R and the last step follows from an expansion near the surface. As it turns out, this is a good approximation 
throughout the star. From this we see that while nr/R is indeed large in the outer regions of the star (the divergence at the 
surface is due to the sound speed vanishing there), it vanishes at the centre and is certainly not large compared to m + 2 in 
fluT) in the inner regions of the star. This is relevant for our analysis since we are trying to work out the induced flow in the 
core, not just the corrections in the thin boundary layer at the base of the crust. 

Given the approximate expression for k it is straightforward to write down an analytic solution to 11011 . If we focus on 
the I = m = 2 r-mode and express the right-hand side as Bx 3 , then the solution is simply 

This should be compared to the incompressible result W-^ nc = Bx 3 /7. Repeating the analysis following 11531 we find that 
s 1 - R c /R 



s 1 - 7R C /8R 



(106) 



This shows that the compressibility lowers the Ekman damping rate by a factor of two or so (for R c /R = 0. 9). Although this 
is far from irrelevant, it is certainly not as dramatic a suppression as the results of lAbnev fc Epstein! <ll99d) suggest. 

The above analysis is, as far as we know, the first attempt to account for the effect of compressibility on the Ekman layer 
damping of r-modes. We have shown how the problem can be solved (up to quadrature), and then estimated the magnitude 
of the effect for a polytropic model. Would it now be useful to consider an even more realistic model based on some tabulated 
equation of state? We think that the answer is no. There are two reasons for this. We have already discussed the well-known 
fact that one should be careful when using realistic (relativistic) equations of state in a Newtonian calculation, and that it is 
difficult to make sure that one is comparing like with like. A second reason why we feel that a more "detailed" calculation may 
not be relevant is that we are still not considering the true physics problem. In a mature neutron star we should incorporate 
both superfluid components at the crust-core interface and the magnetic field. We have shown that the Ekman layer damping 
in an incompressible normal fluid core differs from a compressible model by a factor of a few. This is typical of the uncertainties 
associated with (say) the equation of state, and is perhaps the level of accuracy that we should expect to be able to achieve. 
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7.3 Generalising the surface integral approach 

To complete our discussion of the compressible problem it is relevant to consider how the surface integral approach discussed 
in Section [K] would have to be altered. Writing the Euler equations in terms of enthalpy, as in Section [5] instead of pressure 
we have for the inviscid flow: 

ia v° + 2ft x v° = -V5h° (107) 
Similarly, the 0(^ 1//2 ) induced flow is described by 

icrov 1 + 2fl x v 1 = -VSh 1 + sv° (108) 

We have seen that [s, Sh 1 ] ~ 0(n 1/2 ) and v 1 ~ e>(rT 1/2 ). 
Repeating the steps described in Section |K| we obtain, 

s|v°| 2 = V ■ [v 1 ^ + v'Wt 1 ] - (V ■ v 1 )^ - (V • v )^ 1 (109) 

The first three terms are 0(fi 1/2 ). The last term is much smaller, O{0. 5/2 ), and can therefore be neglected. Transforming the 
first integral to a surface integral in the standard way, we get at leading order, 

Iv ! 2 ^ = /(f ■ v^ShPdS - /(V ■ v^ShPdV (110) 



This is the generalisation of eqn. 1)82^ for the case of a stratified and compressible fluid. The main difference from the 
incompressible result is the presence of the second term in the right-hand side. Since we know that, cf Eq. I|97^ . 

W = -^ (in) 

we see that it encodes the dependence on the compressibility. Since it is a volume integral, it is also clear that the answer 
depends on the compressibility throughout the star not just in the boundary layer. 



8 SUMMARY 

We have revisited the problem of the damping of neutron star r-modes due to the presence of an Ekman layer at the core-crust 
interface. This problem continues to be of interest since there is a possibility that the gravitational- wave driven instability of 
the r-modes is active in rapidly spinning neutron stars, eg. in low-mass X-ray binaries. The Ekman layer is known to provide 
one of the key damping mechanism for these oscillations. 

We reviewed various approaches to the problem and carried out an analytic calculation of the effects due to the Ekman 
layer for a rigid crust. Our obtained estimates support previous numerical results, and provide some useful insights into the 
problem. Of particular interest should be the discussion of the breakdown of the "standard" solution to the problem associated 
with singularities in the velocity field in certain angular directions. Our comparison of previous results, which demonstrates 
that they are all consistent, is also useful. 

We extended the analysis to include the effect of compressibility and composition stratification. This led to a demonstration 
that stratification is unimportant for the r-mode problem, but compressibility suppresses the damping rate by about a factor 
of two (depending on the detailed equation of state). The physical interpretation is that the compressibility counteracts the 
induced flow, and load s to a slightly less efficient energy transport away from the boundary layer. Our result corrects a 
similar analysis of lAbnev fc Epstein! 1^996) for the spin-up problem (associated with pulsar glitches). By assuming that the 
compressibility can be taken to be constant they predicted an exponential suppression of the Ekman layer damping. The 
reason why this assumption is not valid should be clear from our discussion. 

To conclude, we think that this study clears the stage for work on more physically realistic models. The methods required 
to model Ekman layers in rotating neutron stars are now well understood. In particular, we understand the strengths and 
weaknesses of the different routes to the answer. In order to make our models more realistic we need to account for both 
the pres e nce o f ma gnetic fields and superflui d components. The effect of magnetic boundary layers have been discussed by 
iMendelJ l|200l[) and lKinnev fc Mendelll l)2003h . Their work shows that the magnetic field is of key importance. There are also 
many issues associated with superfluid neutron stars, ranging from the presence of additional dynamical degrees of freedom 
to novel dissipation mechanisms. It is clear that, while we have made progress on the simplest Ekman layer problems, many 
challenges remain. 
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APPENDIX 

In this Appendix we provide the steps needed to simplify the angular integral in Section [I] We need to work out 



J = 2n / sin<9Y T 



[ m + l 







im(m + 1) 



1 1 

Ai A2 



YZ-AR c ^-BR c ^}d6 (112) 



To do this we use 



! / I 



sme^Y™Y™ +1 de = - sme~-Y™Y™ +l dO (113) 

Ai J A 2 

and 

sm9A^Y^ +1 de= f S m6B^Y™ +l d6 (114) 

A l JO A 2 

In other words, we have 

J sm ro+1 i — 2x[ — Ym ~ ARc ^fj dd ^ ' 

We can simplify this further in the following way; 

sin 6Y^ +1 AR c ^-d9 = - ( sm6Y™ +1 AR c d e ( i-) dB = [ (sinOY^ARc) d9 

1 / srn^ I + i}Y™ +1 Y™ r^Ym+iY™ - {d e Y™ +1 )(d g Y™) 

2 J Q Ai |_ sm 8 

+ ^MY^+iY™)} dO (116) 
This way we arrive at Eq. I|75^ in Section 0] 
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